function dphidt = phiode(t,phi,eks,P)
dphidt = [(2*pi/(P*(1-eks^2)^(3/2)))*(1+eks*cos(phi(1)))^2];
end